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Abstract 



We report an extensive Monte Carlo study of critical end point behaviour 
in a symmetrical binary fluid mixture. On the basis of general scaling ar- 
guments, singular behaviour is predicted in the diameter of the liquid-gas 
coexistence curve as the critical end point is approached. The simulation re- 
sults show clear evidence for this singularity, as well as confirming a previously 
predicted singularity in the coexistence chemical potential. Both singularities 
should be detectable experimentally. 
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A critical end point occurs when a line of second order phase transitions intersects and is 
truncated by a first order phase boundary delimiting a new noncritical phase. Critical end 
points are common features in the phase diagrams of a diverse range of physical systems, 
notably binary fluid mixtures, superfluids, binary alloys, liquid crystals, certain ferromag- 
nets and ferroelectrics etc. Perhaps the simplest of these is the binary fluid mixture, for 
which the phase diagram is spanned by three thermodynamic fields (T, fi, h), where T is the 
temperature, \l is a chemical potential, and h is an ordering field coupling to the relative 
concentrations of the two fluid components. In the region h = 0, two fluid phases fi and 7 
coexist. By tuning T and fi, however, one finds a critical 'A-line', T c (/x), where both phases 
merge into a single fij phase. This A-line meets the first order line of liquid-gas transitions 
fJL a (T) at a critical end point (T e ,/x e ), see fig 1. For T < T e , the phase boundary fi a (T) 
constitutes a triple line along which the fluid phases fi and 7 coexist with the gas phase a, 
while for T > T e , /i CT (T) defines the region where the and a phases coexist. Precisely 
at the critical end point the critical mixture of (3 and 7 phases coexists with the gas phase. 
Since the gas phase is noncritical, it is commonly referred to as a "spectator" phase. 

Despite their ubiquity, it has only quite recently been pointed out that critical end 
points should exhibit novel properties beyond those observable on the critical line. Using 
phenomenological scaling and thermodynamic arguments, Fisher and coworkers [I] predicted 
that the singular behaviour at the critical end point also engenders new singularities in the 
first order phase boundary itself. Additionally, new universal amplitude ratios were proposed 
for the shape of this boundary, as well as for the noncritical surface tensions near the critical 
end point 0. These predictions were subsequently corroborated by analytical calculations on 
extended spherical models ||. To date, however, empirical support from physically realistic 
systems has been scarce. While some experimental work on surface amplitude ratios has 
been reported [Q, no attention seems to have been given to the bulk coexistence properties 
and the question of the predicted singularity in the spectator phase boundary. There is a 
similar dearth of simulation work on the subject, and we know of no detailed numerical 
studies of critical end point behaviour, either in lattice or continuum models. 

In this letter we move to remedy this situation by providing the first simulation evidence 
for singular behaviour in the first order phase boundary close to a critical end point. The key 
features of our results are as follows. We consider a classical binary fluid within the grand 
canonical ensemble. We review the scaling arguments of Fisher and coworkers and show 
that in addition to the previously predicted singularity in /v(T), they also imply singular 
behaviour in the diameter of the liquid- gas coexistence curve at the critical end point. We 
test these predictions using extensive Monte Carlo simulations of a symmetrical Lennard- 
Jones binary mixture, making full use of modern sampling methods, histogram extrapolation 
techniques and finite-size scaling analyses. The results provide remarkably clear signatures 
of divergences in the appropriate temperature derivatives of the coexistence diameter and 
the phase boundary chemical potential. 

We shall focus attention on the liquid-gas coexistence region in the vicinity of a critical 
end point. Following Fisher and Barbosa M, coexistence is prescribed by the equality of 
the Gibbs free energy G = —IcsTlnE in the gas and liquid phases i.e. G g (fi a (T),T,h) = 
Gi(fj, a (T),T,h). Since the gas spectator phase is noncritical, its free energy is analytic at 
the end point and thus may be expanded as 

G g (fi, T, h) = G e + G{ Afi + G\t + G%h + GIA/j 2 + ■■■ (1) 
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where t = (T — T e )/T e , A/x = /x — /x e < 0. The liquid phase on the other hand is critical 
and therefore contains both an analytic (background) and a singular contribution to the free 
energy 

difi, T, h) = G (/x, T, h) - |r| 2 -^ ± (/i|r|- A ) (2) 

where Go is the analytic part, while Q±{y) is a universal scaling function whose two branches 
± must satisfy matching conditions as y — > ±00. The quantities r(/i, T, h) and h(fi, T, h) are 
both scaling fields that measure deviations from criticality, and comprise linear combinations 
of A/x, t and h. a and A are respectively the specific heat and gap exponents associated 
with the A-line. 

Expanding the critical free energy in A/x, t and h, and setting G g (fi a (T),T, h) = 
Gi{n a (T), T, h), then yields U in the region h — 

fi a (T)-fi (T)^-X ± \t\ 2 - a , T — > T e ± (3) 

with fi (T) = /x e + git H , analytic, X ± > 0. If a > 0, this in turn implies a divergence in 

the curvature of the spectator phase boundary 

^ « -* ± M- (4) 

where the amplitude ratio is expected to be universal J1J. 

Let us now consider the behaviour of the coexistence density in the neighbourhood of 
the critical end point. Specifically, we shall examine the temperature dependence of the 
coexistence diameter in the symmetry plane h — 0, defined by 

Pd{T) = \\pMT)) + pi{iM a {T))) (5) 

which is simply obtained from the coexistence free energy as 

1 f dG(^(T),T) \ 

pdiT) = -v{ dH J (6) 

where G(/x CT (T), T) = [G 5 (/i CT (T), T) + Gi(/i a (T), T)]/2. Appealing to eqs. |l|and||, one then 
finds 

p d (T) = -U ± \r\ 1 - a -V ± \rf 

+ terms analytic at T e (7) 

This singularity is of the same form as the overall density singularity on the critical 
line T c (/x), which for binary fluids with short ranged interactions is expected to be Ising-like. 
For the symmetrical binary fluid studied in the present work, one finds that the amplitude 
V± = in eq. [7[ The diameter derivative then exhibits a specific heat like divergence: 

^ « (8) 
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where we have used |r| = |t|[l + 0(|t| 1_a )] alongthe coexistence curve. Since this divergence 
occurs in the first derivative of Pd{T), it is in principle more readily visible than that in the 
second derivative of p a (T), cf. eq. f|. As we shall now show, however, clear signatures of 
both divergences are readily demonstrable by Monte-Carlo simulation ||. 

The simulations described here were performed for a symmetrical binary fluid model 
using a Metropolis algorithm within the grand canonical ensemble (GCE) 0. The fluid 
is assumed to be contained in (periodic) volume V = L 3 , with grand-canonical partition 
function 

oo oo N ( „ n 

z l = E E II { / e^-^WH^-Ay] (g) 

N 1 =0N 2 =0i=l ^ ' 

where $ = J2i<j 0( r r/) is the total configurational energy, \i is the chemical potential, and h 
is the ordering field (all in units of ksT). N = Ni + N2 is the total number of particles of 
types 1 and 2. The interaction potential between particles i and j was assigned the familiar 
Lennard- Jones (LJ) form 

<P{r lj ) = Ae mn [{a/r l3 ) 12 -{a/r l] )% (10) 

where a is a parameter which serves to set the interaction range, while e mn measures the 
well-depth for interactions between particles of types m and n. In common with most other 
simulations of Lennard- Jones systems, the potential was cutoff at a radius r c = 2.5a to 
reduce the computational effort. 

An Ising model type symmetry was imposed on the model by choosing en = 622 = e > 0. 
This choice endows the system with energetic invariance under h — > —h and ensures that the 
critical end point lies in the symmetry plane h = 0. A further parameter 612 = 8 was used to 
control interactions between unlike particles. The phase diagram of the model in the surface 
h = is then uniquely parameterised by the ratio 8/e. Choosing 8/e < 1 yields a phase 
diagram having a critical end point temperature T e <C T l c 9 and density p e ^> p 1 ^, where T l c 9 
and p l Q are the liquid-gas critical temperature and density respectively. Choosing a smaller 
value of 8/e, however, moves the end point towards the liquid-gas critical point, into which 
it merges for a certain sufficiently small 8/e, forming a tricritical point. Empirically we find 
that the phase diagram is rather sensitive to the choice of 8/e. Thus for 8/e ~ 0.6, we 
find a tricritical point, while for 8/e = 0.75 there is a critical end point having p e w 2.3p^. 
In the present work, all simulations were performed with 8/e = 0.7, which yields critical 
end point parameters T e w 0.93T C ' 9 , p e pa l.Top 1 //. This temperature is sufficiently small 
compared to T l c 9 that critical density fluctuations (which might otherwise obscure the end 
point behaviour) may safely be neglected, while at the same time p e is not so large as to 
hinder particle insertions. 

Although use of conventional GCE simulations to study liquid-gas phase coexistence 
presents no great practical difficulties when T < T l c g [|J, investigations of the strongly 
first order regime are rendered extremely problematic by the large free energy 

barrier separating the coexisting phases. This leads to metastability effects and protracted 
correlation times. To circumvent this difficulty we employed the multicanonical preweighting 
method ||, which encourages the simulation to sample the interfacial configurations of 
intrinsically low probability. This is achieved by incorporating a suitably chosen weight 
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function in the Monte-Carlo update probabilities. The weights are subsequently folded 
out from the sampled distributions to yield the correct Boltzmann distributed quantities. 
Further details of the implementation of this technique as well as a method for determining 
a suitable preweighting function are given elsewhere . 

In the course of the simulations, three systems sizes of volume V = (7.5a) 3 , V = (10cx) 3 
and V = (12.5<r) 3 were studied, corresponding to average particle numbers of N « 250, 
N fa 600 and N « 1200 respectively at the critical end point (whose location we discuss 
below). Following equilibration, runs comprising up to 6 x 10 9 MCS JL0] were performed 
and the density p = N/V, energy density u — <f>/V and number difference order parameter 
m = (Ni — N 2 )/V were sampled approximately every 10 4 MCS. Attention was focused on 
the finite-size distributions Pl(p) and pi^m). Precisely on the liquid-gas coexistence curve, 
the density distribution Pl{p) is (to within corrections exponentially small in L) double 
peaked with equal weight in both peaks [11]. For a given simulation temperature, this 



'equal weight' criterion can be used to determine the coexistence chemical potential to high 
accuracy. Simulations were carried out for each L at several (typically 5) temperatures 
along the coexistence curve, and histogram reweighting [12 1 was used to interpolate between 
simulation points and to aid the precise location of the coexistence chemical potential ||. 
The position of the critical end point itself was estimated using finite-size scaling techniques 
in the standard manner [Tj|, by studying the scaling of the fourth order cumulant ratio 



Ul — 1 — 3(m 4 )/(m 2 ) 2 for pi,(m) as a function of T and L along the liquid branch of the 
coexistence curve. A cumulant intersection || implying critical scale invariance was obtained 
at T e = 0.958(2), p e = 0.587(5), where T = ksTje. Further points on the A-line away from 
the critical end point were also determined using the same method. Related finite-size scaling 
techniques, this time focusing on the density-like ordering operator were utilised to 

locate the liquid-gas critical point. This yielded the estimates T l c 9 = 1.024(2), p l J> = 0.327(2). 

Figure [| shows the estimated coexistence liquid and gas densities as a function of tem- 
perature, determined as the peak densities of Pl{p)- Also shown is the measured locus of 
the A-line and the position of the liquid-gas critical point. Clearly a pronounced 'kink' is 
discernible in the liquid-branch density in the vicinity of the critical end point. The gas 
branch, on the other hand displays no such kink due to the analyticity of G g (p, T) at T e . To 
probe more closely the behaviour of the coexistence density, we plot in figure |3](a) the diam- 
eter derivative —dpd/dT, for the three system sizes studied. The data exhibit a clear peak 
close to T e , that narrows and grows with increasing system size. Very similar behaviour is 
also observed in the curvature of the spectator phase boundary —d 2 p a /dT 2 , see figure 0(b). 
These peaks constitute, we believe, the finite-size-rounded forms of the divergences eqs. |8| 
and [|. On the basis of finite-size scaling theory [13], the peaks are expected to grow in 



height like L a ^ u , with v the correlation length exponent. Unfortunately it is not generally 
feasible to extract estimates of a/v in this way (even for simulations of lattice Ising models), 
because to do so necessitates an accurate measurement of the analytic background, for which 
the present system sizes are much too small. Nevertheless, the correspondence of the peak 
position with the independently estimated value of T e , as well as the narrowing and growth 
of the peak with increasing L constitutes strong evidence supporting the existence of the 
predicted singularities. 

In summary we have employed advanced Monte-Carlo simulation techniques to study the 
first order phase boundary near the critical end point of a continuum binary fluid model. 
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The results provide the first empirical evidence for singularities in the phase boundary and 
the coexistence curve diameter. We expect that similar effects should be experimentally 
observable (not only in binary fluids), and that due to the absence of finite-size rounding 
they should be even more marked than observed here. Moreover, as we have shown, for 
asymmetrical systems such as real binary mixtures, the coexistence diameter is expected to 
manifest a much stronger singularity than occurs in the present symmetrical model. This 
should therefore be particularly conspicuous. 

The author thanks K. Binder and D.P. Landau for stimulating discussions. Helpful 
correspondence with A.D. Bruce, M.E. Fisher, M. Krech and M. Miiller is also gratefully 
acknowledged. This work was supported by BMBF project number 03N8008 C. 
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Wilding figure 1 

FIG. 1. Schematic phase diagram of a binary fluid in the coexistence surface h = 0. The broken 
line Hct(T) is the first order liquid-gas phase boundary terminating at the critical point (cross). 
The full line is the critical line of second order transitions T c (/j.) separating the demixed phases 
(3+j, from the mixed phase (3j. The two lines intersect at the critical end point. 
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Wilding figure 2 

FIG. 2. The peak densities corresponding to the coexistence form of pl(p) for the three systems 
sizes studied, plotted as a function of the temperature. Also shown is the estimated locus of the 
A- line (circles) and the liquid-gas critical point (cross). Statistical errors do not exceed the symbol 
sizes. 
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FIG. 3. (a) The numerical temperature derivative of the measured coexistence diameter 
—dpd/dT* in the vicinity of the critical end point temperature, (b) The measured curvature 
of the phase boundary, —d 2 [i CT /dT 2 , in the vicinity of the critical end point temperature. 
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